function [residual] = kappa_iterative_residual( kp, par ) 
% Usage: kappa_iterative_residual( kp, par ) 
% 
% Note: input kp should be a number. This function returns the residual of the jump equation for \kappa^p.

beta=par.beta;   alpha=par.alpha;   
max_jump = par.max_jump;   xK_roots_fitted=par.xK_roots_fitted;   
w= par.w_ind;      lambda=par.lambda;   
p_2D_fitted=par.p_2D_fitted;     % p_2D_fitted: p_fitted is an object generated by scatteredInterpolant

p = p_2D_fitted(  w, lambda );
xK = ppval( xK_roots_fitted, kp);  
delta_x = beta*max(xK-1,0);
kappab= (xK*kp+alpha*delta_x);
kappa_fs = alpha*delta_x * w / (1-w);
yK = (1-w*xK)/(1-w);
kappah = par.varepsilon + (1-par.varepsilon)*yK*kp - kappa_fs;

w_new = w .* (  1- kappab ) ./  ( 1- w*kappab - (1-w)*kappah )  ;
lambda_new = lambda + par.kappa_lambda_fun(lambda);

residual = 1 - p_2D_fitted( w_new, lambda_new)/p - kp  + 10^4*(kp<0 | kp>max_jump) ;
